{
  "metadata": {
    "kernelspec": {
      "name": "xpython",
      "display_name": "Python 3.13 (XPython)",
      "language": "python"
    },
    "language_info": {
      "name": ""
    }
  },
  "nbformat_minor": 5,
  "nbformat": 4,
  "cells": [
    {
      "id": "ca6d609a-fb57-4d8a-97a7-ae2c4a09baa1",
      "cell_type": "code",
      "source": "import numpy as np\nimport pandas as pd\nfrom qutip import Qobj, sigmax, sigmay, sigmaz, mesolve, Options\nfrom scipy.optimize import curve_fit\nimport statsmodels.api as sm\n\n# Define exponential decay model\ndef exp_decay(t, a, tau, c):\n    return a * np.exp(-t / tau) + c\n\n# Simulation parameters\nomegas = np.array([0.05e6, 0.1e6, 0.2e6, 0.5e6])  # Rabi frequencies (MHz)\ngammas = np.array([0.05e6, 0.1e6, 0.2e6])          # Decoherence rates (MHz)\ntlist = np.linspace(0, 50e-6, 500)  # Time array (µs)\n\n# TLS Hamiltonian components\nsigma_x, sigma_y, sigma_z = sigmax(), sigmay(), sigmaz()\nH0 = 0.5 * 2 * np.pi * 3e8 / 630e-9 * sigma_z  # Transition frequency from 630 nm wavelength\nH1 = 0.5 * sigmax()\n\n# Function to simulate and fit coherence decay\ndef fit_coherence_decay(omega, gamma):\n    H = H0 + 0.5 * omega * H1\n    psi0 = Qobj([[1], [0]])\n    c_ops = [np.sqrt(g) * s for g, s in zip([gamma], [sigma_x, sigma_y, sigma_z])]\n    result = mesolve(H, psi0, tlist, c_ops=c_ops, e_ops=[sigmax(), sigmay()], options=Options(store_states=True))\n    coherence = np.sqrt(result.expect[0]**2 + result.expect[1]**2)\n    \n    # Fit exponential decay\n    popt, pcov = curve_fit(exp_decay, tlist, coherence, p0=[1.0, 1e-6, 0.0], maxfev=10000)\n    a, tau, c = popt\n    \n    # Statistical analysis for p-values (using OLS on residuals)\n    residuals = coherence - exp_decay(tlist, *popt)\n    X = np.vstack([np.ones(len(tlist)), np.exp(-tlist / tau), tlist * 0 + 1]).T  # Model matrix\n    model = sm.OLS(residuals, X).fit()\n    p_values = model.pvalues\n    \n    return [omega * 1e-6, gamma * 1e-6, a, tau * 1e6, c, p_values[1], p_values[2], p_values[0]]  # τ in µs\n\n# Generate fit dataset\nfit_data = []\nfor omega in omegas:\n    for gamma in gammas:\n        fit_params = fit_coherence_decay(omega, gamma)\n        fit_data.append(fit_params)\n\n# Save to CSV\nfit_df = pd.DataFrame(fit_data, columns=['Omega (MHz)', 'Gamma (MHz)', 'Amplitude (a)', 'Decay_Time (τ µs)', 'Offset (c)', 'P_Value_a', 'P_Value_τ', 'P_Value_c'])\nfit_df.to_csv('tls_coherence_decay_fits.csv', index=False)",
      "metadata": {
        "trusted": true
      },
      "outputs": [],
      "execution_count": null
    },
    {
      "id": "a1b66e05-74d3-4959-befc-8717084cf17e",
      "cell_type": "code",
      "source": "",
      "metadata": {
        "trusted": true
      },
      "outputs": [],
      "execution_count": null
    }
  ]
}